Magnetic Quenching of Turbulent DifFusivity: Reconciling 
Mixing- length Theory Estimates with Kinematic Dynamo Models 

of the Solar Cycle 

Andres Mufioz-Jaramillo 

Department of Physics, Montana State University, Bozeman, MT 59717, USA 

O munozOsolar .physics .montana.edu 

(N 

'^ and 

I^ Dibyendu Nandy 

n/! Indian Institute for Science Education and Research, Kolkata, Mohampur 741252, West 

C/^ Bengal, India 

Ph dnandi@iiserkol.ac. in 

6 



c/5 



o 
o 

> 

X 

S3 



and 
Petrus C. H. Martens 



J> Department of Physics, Montana State University, Bozeman, MT 59717, USA 

(N 
^ Harvard- Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA 



martensOsolar . physics . montana . edu 

ABSTRACT 

The turbulent magnetic diffusivity in the solar convection zone is one of the 
most poorly constrained ingredients of mean-field dynamo models. This lack 
of constraint has previously led to controversy regarding the most appropriate 
set of parameters, as different assumptions on the value of turbulent diffusivity 
lead to radically different solar cycle predictions. Typically, the dynamo commu- 
nity uses double step diffusivity profiles characterized by low values of diffusivity 
in the bulk of the convection zone. However, these low diffusivity values are 
not consistent with theoretical estimates based on mixing-length theory - which 
suggest much higher values for turbulent diffusivity. To make matters worse, 
kinematic dynamo simulations cannot yield sustainable magnetic cycles using 



these theoretical estimates. In this work we show that magnetic cycles become 
viable if we combine the theoretically estimated diffusivity profile with magnetic 
quenching of the diffusivity. Furthermore, we find that the main features of this 
solution can be reproduced by a dynamo simulation using a prescribed (kine- 
matic) diffusivity profile that is based on the spatiotemporal geometric-average 
of the dynamically quenched diffusivity. Here, we provide an analytic fit to the 
dynamically quenched diffusivity profile, which can be used in kinematic dynamo 
simulations. Having successfully reconciled the mixing-length theory estimated 
diffusivity profile with kinematic dynamo models, we argue that they remain a 
viable tool for understanding the solar magnetic cycle. 

Subject headings: Sun: dynamo. Sun: interior. Sun: activity 



1. Introduction 

The solar magnetic cycle involves the recycling of the toroidal and poloidal components 
of the magnetic field which are generated at spatially segregated source layers that must 
communicate with each-other (see e.g., Wilmot-Smith et al. 2006; Charbonneau 2005). This 
communication is mediated via magnetic flux-transport, which in most kinematic solar dy- 
namo models, is achieved through diffusive and advective (i.e., by meridional circulation) 
transport of magnetic fields. The relative strength of turbulent diffusion and meridional 
circulation determines the regime in which the solar cycle operates, and this has far reaching 
implications for cycle memory and solar cycle predictions (Yeates, Nandy & Mackay 2008; 
Nandy 2010). As shown in Yeates, Nandy & Mackay (2008), different assumptions on the 
strength of turbulent diffusivity in the bulk of the Solar Convection Zone (SCZ) lead to differ- 
ent predictions of the solar cycle (Dikpati, DeToma & Oilman 2006; Choudhuri, Chatterjee 
& Jiang 2007). Previously this lack of constraint has led to controversy regarding what value 
of turbulent diffusivity is more appropriate and yields better solar like solutions (Nandy & 
Choudhuri 2002, Dikpati et al. 2002, Chatterjee et al. 2004, Dikpati et al. 2005, Choudhuri 
et al. 2005). Currently, most dynamo modelers use double-step diffusivity profiles which 
are somewhat ad- hoc and different from one-another (see Figure 1; Rempel 2006, Dikpati 
and Gilman 2007, Guerrero and de Gouveia Dal Pino 2007, Jouve and Brun 2007). There 
is however, a way of theoretically estimating the radial dependence of magnetic diffusivity 
based on Mixing Length Theory (MLT; Prandtl 1925). 



2. Order of Magnitude Estimation 

Going back to the derivation of tlie mean-field dynamo equations (after using tlie first 
order smoothing approximation), we find that the turbulent diffusivity coefficient becomes 
(Moffat 1978): 

ri=l{v'), (1) 

where r is the eddy correlation time and v corresponds to the turbulent velocity field. In order 
to make an order of magnitude estimation we turn to MLT, which although not perfect, has 
been found to be in general agreement with numerical simulations of turbulent convection 
(Chan & Sofia 1987; Abbett et al. 1997). More specifically we use the Solar Model S 
(Chistensen-Dalsgaard et al. 1996), which is a comprehensive solar interior model used by 
GONG in all their helioseismic calculations. Among other quantities, this model estimates 
the mixing length parameter a^, the convective velocity v for different radii and the necessary 
variables to calculate the pressure scale height Hp. In terms of those quantities the diffusivity 
becomes: 

V ~ T^apHpV, (2) 

which we plot in Figure ^ (solid black line) and show how it compares to commonly used 
diffusivity profiles. 



3. The Problem and a Possible Solution 

It is evident that there is a major discrepancy between the theoretical estimate and the 
typical values used inside the convection zone (around two orders of magnitude difference), 
dynamo models simply cannot operate under such conditions. 

A possible solution to this inconsistency resides in the back-reaction that strong mag- 
netic fields have on velocity fields, which results in a suppression of turbulence and thus of 
turbulent magnetic diffusivity (Roberts & Soward 1975). This magnetic "quenching" of the 
turbulent diffusivity has been studied before in different contexts (Riidiger et al. 1994; Tobias 
1996; Oilman & Rempel 2005; Muiioz-Jaramillo, Nandy & Martens 2008; Guerrero, Dikpati 
& de Gouveia Dal Pino 2009). However, although this issue has been common knowledge 
for more than a decade, it's only because of current improvements in computational tech- 
niques (Hochbruck & Lubich 1997; Hochbruck, Lubich & Selhofer; Munoz-Jaramillo, Nandy 
& Martens 2009; MNM09 from here on), that this question can be finally addressed quanti- 
tatively. In this paper we study whether introducing magnetic quenching of the diffusivity 
can solve this discrepancy and whether the shape of the currently used diffusivity profiles 



can be understood as a spatiotemporal average of the effective turbulent diffusivity after 
taking quencliing into account. 



4. The Kinematic Mean-Field Dynamo Model 

Our model is based one the axisymmetric dynamo equations: 

(3) 
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where A is the 0-component of the potential vector (from which B.^. and Bq can be obtained), 
B is the toroidal field {B^)^ Vp is the meridional fiow, Q the differential rotation, r] the 
turbulent magnetic diffusivity and s = r sin(^). In order to integrate these equations we need 
to prescribe four ingredients: meridional flow, differential rotation, poloidal fleld regeneration 
mechanism and turbulent magnetic diffusivity. We use the same meridional flow proflle we 
deflned in MNM09, which better captures the features present in helioseismic data. Our flow 
proflle has a penetration depth of .675Rq and a top speed of 12m/s. For the differential 
rotation we use the analytical form of Charbonneau et al. (1999), with a tachocline centered 
at 0.7Rq and whose thickness is O.O5i?0. For the poloidal fleld regeneration mechanism 
we use the improved ring-duplet algorithm described in Nandy, Muhoz-Jarmillo & Martens 
2010 (NMMIO) and Muhoz-Jaramillo et al. 2010 (MNMYIO), but using a value of Kq = 
3900, in order to insure super-criticality. As in NMMIO and MNMYIO, we restrict active 
region emergence to latitudes between 45°N and 45° S. Speciflcs about our treatment of the 
turbulent diffusivity are described below. More details regarding kinematic dynamo models 
can be found in a review by Charbonneau (2005) and references therein. 



4.1. Turbulent Magnetic Diffusivity and Diffusivity Quenching 

In order to study the effect of magnetic quenching on dynamo models we introduce an 
additional state variable rjmq governed by the following differential equation: 

2/ n .MU2 -Vn,g{r,9,t)\ . (4) 



8t T \1 + B^{r,9,t)/B, 



In a steady state, rjmq corresponds to the MLT estimated diffusivity rjMLrij') quenched in such 
a way that the diffusivity is halved for a magnetic field of amplitude Bq = 6700 Gauss (G). 
This value corresponds to the average equipartition field strength inside the SCZ calculated 
using the Solar Model S. The characteristic time of relaxation r = 30 days is an estimate of 
the average eddy turnover time. 

We make a fit of riMLT{f) using the following analytical profile (see Figure |2]): 

where tjq = lO^cm^/s corresponds to the diffusivity at the bottom of the computational 
domain; rji = 1.4 x lO^^cm'^/s and r/2 = lO^^crn^/s control the diffusivity in the convection 
zone; ri = O.71-R0, di = 0.015Rq, r2 = 0.96Rq and d2 = O.O9i?0 characterize the transitions 
from one value of diffusivity to the other. With this in mind, we define the effective diffusivity 
at any given point as 

rjeffir, 9, t) = r]min{r) + r]mq{r, 9, t). (6) 

with the minimum magnetic diffusivity ?7mm('") given by the following analytical profile (see 
Figure g): 

Vn.Ur) =Vo + ^^^ (l + erf (^) ) , (7) 

where rjcz = 10^°cm^/s, r^z = 0.69Rq, and dcz = O.O7i?0. Since diffusivity is now a state 
variable, small errors can lead to negative values of diffusivity, which in turn leads to unbound 
magnetic field growth. By putting a limit on how small can the diffusivity become, we 
successfully avoid this type of computational instability. 



5. Dynamo Simulations 

We use the SD-Exp4 code (see the Appendix of MNM09) to solve the dynamo equations 
(Eqs. Is]). Our computational domain comprises the SCZ and upper layer of the solar radiative 
zone in the northern hemisphere (O.SS-Rq < r < Rq and < 9 < 7c). In order to approximate 
the spatial differential operators with finite differences we use a uniform grid (in radius and 
co-latitude), with a resolution of 400 x 400 gridpoints. 

Our boundary conditions assume that the magnetic field is anti-symmetric across the 
equator (dA/d9\d=T,/2 = 0; dB /d9\e=n/2 = 0), that the plasma below the lower boundary is 
a perfect conductor {A{r = 0.55Rq,9) = 0; d{rB) / dr\r=o.55RQ = 0), that the magnetic field 
is axisymmetric {A{r,9 = 0) = 0; B{r,9 = 0) = 0), and that field at the surface is radial 
{d{rA)/dr\r=RQ = 0; B{r = Rq,9) = 0). Our initial conditions consist of a large toroidal 



belt and no poloidal component. After setting up the problem we let the magnetic field 
evolve for 200 years allowing the dynamo to reach a stable cycle. 



6. Results and Discussion 

The first important result is the existence of a uniform cycle in dynamic equilibrium. 
The presence of a diffusivity quenching algorithm allows the dynamo to become viable in a 
regime in which kinematic dynamo models cannot operate thanks to the creation of pockets 
of relatively low magnetic diffusivity (where long lived magnetic structures can exist). This 
can be clearly seen in Figure [3| which shows snapshots of the effective turbulent diffusivity 
and the toroidal and poloidal components of the magnetic field at different moments during 
the sunspot cycle (half a magnetic cycle). As expected the turbulent diffusivity is strongly 
suppressed by the magnetic field (especially by the toroidal component), increasing the 
diffusive timescale to the point where diffusion and advection become equally important for 
fiux transport dynamics. This slow-down of the diffusive process is crucial for the survival of 
the magnetic cycle since it gives differential rotation more time to amplify the weak poloidal 
components of the magnetic field into strong toroidal belts, while providing them a measure 
of isolation from the top (r = Rq) and polar {6 = 0) boundary conditions (B=0). 



6.1. Time and Spatiotemporal Averages of the Effective Diffusivity 

Given that ultimately we want to understand how adequate kinematic diffusivity profiles 
are and whether they are plausible representations of physical reality consistent with MLT, 
we need to find a connection between kinematic profiles and the dynamically quenched 
diffusivity. Because of this, the next natural step is to find time and spatiotemporal averages 
of the effective diffusivity. After careful consideration we have found that the geometric 
average (also known as logarithmic average): 

Nt 

^og[r]avT{ri,Oj)] = j^^log[r/ej/(ri,ej,t„)] 

71=0 

(8) 

Ng 

j=0 

shown in Figure Sa (time average) and as a solid line in Figure Sc (spatiotemporal aver- 
age), captures better the physical and mathematical nuances of diffusive processes than the 



arithmetic average: 

Nt 
n=0 

(9) 

Ne 
VavTL{ri) = ^^^VavT{ri,ej), 
j=0 

shown in Figure Sb (time average) and as a sohd hne in Figure Sc (spatiotemporal average). 
The reason is that the geometric average is more appropriate for systems whose evolution 
is exponential in nature and when the values involved in the average are dependent of each 
other. Since the diffusion equation has exponential solutions and our intention is to capture 
the essence of a cyclic behavior in which the magnetic field is diffused and advected on a 
closed circuit (thus there is a direct dependence between the different elements considered 
in our average), the geometric average is the best choice. As an additional advantage, the 
geometric mean gives more weight to small values than the arithmetic mean; this is important 
because by definition the regions of depressed diffusivity are also those which the magnetic 
field visits during the cycle. Figure [3] shows that turbulent diffusivity is being suppressed 
mainly in a region centered at mid-latitudes, whereas the polar and equatorial regions remain 
unquenched (and devoid of magnetic field). If one then compares the geometric time average 
(Fig. El-a) and the arithmetic time average (Fig. El-a) with the evolution of the magnetic 
diffusivity (Fig. [3]), it's clear that the geometric average is qualitatively truer to the essence 
of the process. As we will see below, the quantitative results are very encouraging as well. 



6.2. Comparison with Kinematic Dynamo Simulations 

Once we calculate the time and latitude geometric average of the effective diffusivity 
we obtain a radial diffusivity profile (see solid line on Fig. |4]-c), which interestingly can be 
accurately described as a double-step profile (Eq. |5]) with the following parameters: r/o = 
lO^cmVs, 771 = 1.6 X lO^cmVs, 772 = 3.25 x lO^^cm'^/s, ri = 0.71Rq, di = 0M7Rq, 
r2 = O.895-R0 and d2 = O.OSli?© (see circles on Fig. |4]-c). We then run a kinematic dynamo 
simulation (no magnetic quenching), leaving all ingredients intact with the exception of the 
turbulent diffusivity profile - for which we use this double-step fit. 

In order to compare the general properties of both simulations we cast the results in 
the shape of synoptic maps (also known as butterfly diagrams) as can be seen in Figure [s] 
The results obtained using the MLT estimate and diffusivity quenching (Fig. (SJ-a) and the 
results obtained using a kinematic simulation with the geometric average fit (Fig. [SJ-c) are 
remarkably similar given the very different nature of the two simulations. It is clear that the 



shape of the solutions differs mainly in the active region emergence pattern. However, the 
general properties of the cycle (amplitude, period and phase) are successfully captured by the 
geometric average and are essentially the same. This result argues in favor of the capacity 
of kinematic diffusivity profiles of capturing the essence of turbulent magnetic quenching. 
These results define a framework which can be used to find a physically meaningful diffusivity 
profile based upon fundamental theory and models of the solar interior, rather than through 
heuristic approaches. 



6.3. Comparison with Solar Cycle Observations 

Ultimately, the goal of dynamo models is to understand the solar magnetic cycle and 
reproduce and predict it's main characteristics. It's therefore important to compare our 
results with solar observations. It is clear that the solutions are not exactly similar to those 
of kinematic dynamo simulations whose parameters have been finely tuned: cycle period 
of 7 years instead of 11, broad wings and incorrect phase. This differences point to an 
overestimation of the turbulent diffusivity; mainly near the surface (affecting phase and 
period) and at the bottom of the SCZ (which affects period and the shape of the wings). 
The cause of this overestimation resides in our definition of diffusivity quenching: in this 
work we use the average kinetic energy present in convection, which means that diffusivity 
is quenched equally through the convection zone. However, convection is less energetic near 
the bottom of the SCZ (due to low convective speeds) and near the surface (due to low 
mass density). This means that simulations taking this factor into account will probably 
yield more correct solutions. Nevertheless, the solutions we obtain are reasonably accurate, 
given the fact that we have completely refrained from doing any fine tuning, instead limiting 
ourselves to fundamental theories and models. 



7. Conclusions 

In summary, we have shown that coupling magnetic quenching of turbulent diffusivity 
with the estimated profile from mixing length theory, allows kinematic dynamo simulations 
to produce solar-like magnetic cycles, which was not achieved before. Therefore, we have 
reconciled mixing length theory estimates of turbulent diffusivity with kinematic dynamo 
models of the solar cycle. Additionally, we have demonstrated that kinematic simulations 
using a prescribed diffusivity profile based on the geometric average of the dynamically 
quenched turbulent diffusivity, are able to reproduce the most important cycle characteristics 
(amplitude, period and phase) of the non-kinematic simulations. Incidentally, this radial 



profile can be described by a double step profile, which has been used extensively in recent 
solar dynamo simulations. From the simulations reported here we provide an analytic fit to 
this double-step diffusivity profile that best captures the effect of magnetic quenching. A 
posteriori, our results strongly support the use of kinematic dynamo simulations as tools for 
exploring the origin and variability of solar magnetic cycles. 
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Radial Dependence of the Turbulent Magnetic Diffusivity 
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Fig. 1. — Different diffusivity profiles used in kinematic dynamo simulations. The solid 
black line corresponds to an estimate of turbulent diffusivity obtained by combining Mixing 
Length Theory (MLT) and the Solar Model S. The fact that viable solutions can be obtained 
with such a varied array of profiles have led to debates regarding which profile is more 
appropriate. Nevertheless, it is well known that kinematic dynamo simulations cannot yield 
viable solutions using the MLT estimate. 
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Fig. 2. — Fit (solid line) of diffusivity as a function of radius to the mixing- length theory 
estimate (circles). As part of our definition of effective diffusivity we put a limit on how 
much can the diffusivity be quenched. This minimum diffusivity has a radial dependence 
that is shown as a dashed line. 
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Fig. 3. — Snapshots of the effective diffusivity and the magnetic field over half a dynamo 
cycle (a sunspot cycle). For the poloidal field a solid (dashed) line corresponds to clockwise 
(counter-clockwise) poloidal field lines. Each row is advanced in time by an sixth of the 
dynamo cycle (a third of the sunspot cycle) i.e., from top to bottom t = 0, r/6, r/3 and r/2. 
As expected, the turbulent diffusivity is strongly depressed by the magnetic field (especially 
by the toroidal component). This reduces the diffusive time-scale to a point where the 
magnetic cycle becomes viable and sustainable. 
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Fig. 4. — Geometric (Eq. pi) and arithmetic (Eq. ^ averages of the effective diffusivity. We 
find that the geometric time average (a) and geometric spatiotemporal average (sohd hne 
in c) capture the essence of diffusivity quenching much more accurately than the arithmetic 
time average (b) and arithmetic spatiotemporal average (dashed hne in c). We then use a 
fit to the geometric spatiotemporal average (circles in c), in a kinematic dynamo simulation 
in order to find out how it compares with the simulation that uses the MLT estimate and 
magnetic quenching. 
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Surface Radial Field and AR emergence 
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Fig. 5. — Synoptic maps (butterfly diagrams) showing the time evolution of the magnetic fleld 
in a simulation using the Mixing-Length Theory (MLT) estimate and diffusivity quenching 
(a), and a kinematic simulation using the geometric spatiotemporal average of the dynami- 
cally quenched diffusivity (b). They are obtained by combining the surface radial field and 
active region emergence pattern. For diffuse color, red (blue) corresponds to positive (neg- 
ative) radial field at the surface. The each red (blue) dot corresponds to an active region 
emergence whose leading polarity has positive (negative) fiux. We can see that a kinematic 
dynamo simulation in which we leave all parameters the same, but fix the diffusivity to the 
geometric spatiotemporal average, can capture the most important features of the magnetic 
cycle produced by the simulation using the MLT estimate and diffusivity quenching (period, 
amplitude and phase). 



